#ifndef VORONOI_FUNCTION
#define VORONOI_FUNCTION


extern "C" bool cuPrintInit();
extern "C" void initCUDA(int argc,char** argv);
extern "C" void voronoiClean();
__global__ void vecAdd(float* A , float * B , float *C);
__global__ void voronoi_kernel(float2* d_sitesOrigin ,int* d_sitesOriginLength,Event* d_eventHeap,int* d_eventHeapSize,float2* d_vertexPool,int* d_vertexPoolLength,float2* d_center,ArcHeapEle* d_arcHeap,int* d_arcHeapSize,EllipseArc* d_arcsPool , int* d_arcsPoolLength,ActiveEdge* d_activeEdgePool,int* d_activeEdgePoolLength,int2* d_edgePool,int* d_edgePoolLength,int*  iterateTimes,int* d_headArcsPos,EllipseArcNoPt* d_arcsNoPtPool,int* d_arcsNoPtPoolLength,ActiveEdgeNoPt* d_activeEdgeNoPtPool,int* d_activeEdgeNoPtPoolLength,EventNoPt* d_eventHeapNoPt ,int* d_eventHeapNoPtSize );
__device__ void removeOutofdateEvents(Heap<Event>& eventHeap);
__device__ Event GetNextEvent(Heap<Event>& eventHeap);
void copyArc();
void copyActiveEdge();
void voronoi_kernel_c();

extern "C" bool cuPrintInit()  
{  
	cudaError_t err = cudaPrintfInit();  
	if(0 != strcmp("no error", cudaGetErrorString(err))) {
		puts("device is not ready!\n");
		exit(EXIT_FAILURE);
	}
	return true;
}



//void 
extern "C" void initCUDA(int argc,char** argv){

	
	//init device		
	cudaDeviceProp deviceProp;
	cudaGetDeviceProperties(&deviceProp,0);
	if(!deviceProp.canMapHostMemory){
		fprintf(stderr,"Device %d cannot map host memory!\n",0);
		exit(EXIT_FAILURE);
	}
	cudaSetDeviceFlags(cudaDeviceMapHost);
	

	int pointNum;

	int2 blockGrid;
	int2 threadGrid;
	
	radiusVaring = 0.0f;

	printf("please input  point num(max size %d) :",EVENT_HEAP_SIZE );

#ifdef _SCANF_
	scanf( "%d" , &pointNum);
#else 
	#ifdef _COMMAND_LINE_
		if( argc < argcNum  ) {
			printf("too few argument\n");
			exit(1);
		}
		sscanf(argv[1],"%d",&pointNum);
	#else
		pointNum = 9000;
	#endif
#endif

	do{
		
		printf("please input block grid x and thread grid x (threadgrid y == x):");
#ifdef _SCANF_
		scanf( "%d%d" , &blockGrid.x , &threadGrid.x);
#else
#ifdef _COMMAND_LINE_
		if( argc < argcNum ){
			printf("too few argument\n");
			exit(1);
		}
		sscanf(argv[2],"%d",&blockGrid.x);
		sscanf(argv[3],"%d",&threadGrid.x);
	#else
		blockGrid.x = 8;
		threadGrid.x = 6;
	#endif
#endif
		blockGrid.y = blockGrid.x;
		threadGrid.y = threadGrid.x;
		blockNum = blockGrid.x * blockGrid.y;
		threadNum = threadGrid.x * threadGrid.y;
	}while( threadNum < 1);



	//int sitesMemsize = (((sitesInt.size()*sizeof(float2)+PAGE_SIZE-1)/PAGE_SIZE)*PAGE_SIZE);
	//int sitesLengthMemsize = (((sizeof(int)+PAGE_SIZE-1)/PAGE_SIZE)*PAGE_SIZE);
	//printf("sitesMemsize siteslengthmemsize %d %d\n",sitesMemsize,sitesLengthMemsize);

	//h_sites=(float2*)VirtualAlloc(NULL,sitesMemsize,MEM_COMMIT,PAGE_READWRITE);
	//h_sitesLength=(int*)VirtualAlloc(NULL,sitesLengthMemsize,MEM_COMMIT,PAGE_READWRITE);

	sites.initial(SITES_SIZE,threadNum*blockNum);
	vertices.initial(VERTEX_POOL_SIZE,threadNum*blockNum);
	edges.initial(EDGE_POOL_SIZE,threadNum*blockNum);
	arcs.initial(ARC_POOL_SIZE,threadNum*blockNum);
	arcsNoPt.initial(ARC_POOL_SIZE,threadNum*blockNum);
	centers.initial(1,threadNum*blockNum);
	distances.initial(1,threadNum*blockNum);
	headArcsPos.initial(1,threadNum*blockNum);
	for(int i = 0; i < threadNum*blockNum; ++i)
		headArcsPos.atHost(i,0) = -1;

	bounds.initial(2,threadNum*blockNum);//leftdown and upright
	iterateTimes.initial(1,threadNum*blockNum);
	activeEdges.initial(ACTIVE_EDGE_POOL_SIZE,threadNum*blockNum);
	activeEdgesNoPt.initial(ACTIVE_EDGE_POOL_SIZE,threadNum*blockNum);
	eventHeaps.initial(EVENT_HEAP_SIZE,threadNum*blockNum);
	eventHeapsNoPt.initial(EVENT_HEAP_SIZE,threadNum*blockNum);
	
	arcHeaps.initial(ARC_HEAP_SIZE,threadNum*blockNum);

	struct Point{
		int x,y;
		Point(){}
		Point(float _x,float _y):x(_x),y(_y){}
		bool operator<(const Point& other) const{
			if( x != other.x )
				return x < other.x;
			else 
				return y < other.y;
		}
	};

	float radiusRate;

#ifdef _SCANF_
	printf("please input radiusRate  and sleeptime");
	scanf("%f%d" , &radiusRate,&sleepTime);
	printf("\n");
#else
	#ifdef _COMMAND_LINE_
		if( argc < argcNum ){
			printf("too few argument\n");
			exit(1);
		}
		sscanf(argv[4],"%f",&radiusRate);
		sscanf(argv[5],"%d",&sleepTime);
	#else
		radiusRate = 1.5;
		sleepTime = 10;
	#endif
#endif

	//float tmpDisFromCenterScale;
	int precisionNum;

#ifdef _SCANF_

	//printf("input the scale");
	//scanf("%f" , &tmpDisFromCenterScale );
	printf("input precision num");
	scanf("%d" , &precisionNum );
#else
	#ifdef _COMMAND_LINE_
		if( argc < argcNum ){
			printf("too few argument\n");
			exit(1);
		}
		sscanf(argv[6],"%d",&precisionNum);
	#else
		//tmpDisFromCenterScale = 500;
		precisionNum = 1000;
	#endif

#endif


	//for(int i = 0; i <blockNum * threadNum; ++i){
	//	iterateTimes.atHost(i,0) = tmpIterateTimes;
	//}

	std::tr1::mt19937 eng;
	std::tr1::uniform_int<int> dist(0,precisionNum);
	std::set<Point>rawSites;
	float tmpMinDis = FLT_MAX;

	while(1){
		Point tmp(dist(eng),dist(eng));
		bool filtFlag = false;
#ifdef _FILTE_NEAR_POINT_
		for(std::set<Point>::iterator itr = rawSites.begin(); itr != rawSites.end();++itr){
			if( calcDistance(make_float2(tmp.x,tmp.y),make_float2(itr->x,itr->y)) < WXN_EPSILON * tmpDisFromCenterScale ){
				printf("filter !\n");
				printf("tmpmindis scale %f\n" , tmpMinDis / WXN_EPSILON );
				filtFlag = true;

			}
			tmpMinDis = min(tmpMinDis,calcDistance(make_float2(tmp.x,tmp.y),make_float2(itr->x,itr->y)));
		}
#endif
		if( filtFlag ) continue;
		rawSites.insert(tmp);
		printf("{%f,%f},",(float)tmp.x / (float)precisionNum - .5f,(float)tmp.y/(float)precisionNum - .5f);
		if( rawSites.size() >= pointNum ) break;
	}



	float2 gridLength = make_float2(totalLength/blockGrid.x/(float)threadGrid.x,totalLength/blockGrid.y/(float)threadGrid.y);
		
	for(int bi = 0; bi < blockGrid.x;++bi){
		for(int bj = 0; bj < blockGrid.y;++bj){
			printf("bi %d bj %d\n" , bi , bj );
			for(int i = 0; i < threadGrid.x; ++i){
				for( int j = 0; j < threadGrid.y; ++j){
					int blockId = bi * blockGrid.y + bj;
					int threadId = blockId * threadGrid.x * threadGrid.y +	i*threadGrid.y+j;
					int curI = bi * threadGrid.x + i;
					int curJ = bj * threadGrid.y + j;


					bounds.atHost(threadId,0) = make_float2( gridLength.x * (float)curI - totalHalf, gridLength.y * (float)curJ - totalHalf);
					bounds.atHost(threadId,1) = make_float2( gridLength.x * (float)(curI+1)-totalHalf , gridLength.y * (float)(curJ+1)-totalHalf );
					//printf(" threadId %d bounds %f %f %f %f\n" , threadId, bounds.atHost(threadId,0).x,bounds.atHost(threadId,0).y,bounds.atHost(threadId,1).x,bounds.atHost(threadId,1).y  );
					centers.atHost(threadId,0) = make_float2( ( bounds.atHost(threadId,0).x + bounds.atHost(threadId,1).x ) / 2.0f , (bounds.atHost(threadId,0).y + bounds.atHost(threadId,1).y ) / 2.0f );

					for(std::set<Point>::iterator itr = rawSites.begin(); itr != rawSites.end();++itr ){
						float2 pt = make_float2((float)itr->x / (float)precisionNum - .5f,(float)itr->y / (float)precisionNum - .5f );
						if( calcDistance(pt,centers.atHost(threadId,0)) <= gridLength.x / sqrt(2.0f) * radiusRate ){
							sites.insert(threadId,pt);
							printf( "threadid %d initial Point pt.x %f pt.y %f dis %f scale %f grid banjing %f\n",threadId,pt.x,pt.y,	calcDistance(pt,centers.atHost(threadId,0)),WXN_EPSILON * tmpDisFromCenterScale,gridLength.x / sqrt(2.0f) * radiusRate);
							//std::cout << pt.x << " " << pt.y << "\n";
						}
					}
					if( blockId == 0 && i*threadGrid.y+j == 5 ){
						printf("bounds 0 %f %f\n" , bounds.atHost(threadId,0).x , bounds.atHost(threadId,0).y);
						printf("bounds 1 %f %f\n" , bounds.atHost(threadId,1).x , bounds.atHost(threadId,1).y);
						printf("num %d\n" , sites.length(threadId) );
						printf( "%f %f",centers.atHost(threadId,0).x,centers.atHost(threadId,0).y );
						for(int i = 0; i < sites.length(threadId);++i){
							printf( ",{%f,%f}",sites.atHost(threadId,i).x,sites.atHost(threadId,i).y );
						}
						sites.length(threadId) = 0;
					}

				}
			}
		}
	}

		Event evt;
		for(int threadId = 0;threadId < blockNum * threadNum;++threadId){
			Heap <Event>tmpEventHeap;
			tmpEventHeap.initial(&((eventHeaps.hostPoint())[threadId*EVENT_HEAP_SIZE]),(eventHeaps.hostLength())[threadId],EVENT_HEAP_SIZE);
			for(int i = 0; i < sites.length(threadId);++i){
				evt.type = NewEllipse;
				float2 tmpSite = sites.atHost(threadId,i);
				float2 tmpCenter = centers.atHost(threadId,0);
				evt.site = tmpSite;
				evt.distance = sqrt((tmpCenter.x-tmpSite.x)*(tmpCenter.x-tmpSite.x)+(tmpCenter.y-tmpSite.y)*(tmpCenter.y-tmpSite.y));		
				evt.leftActiveEdge = evt.rightActiveEdge = NULL;
				evt.arc = NULL;
				tmpEventHeap.push(evt);
			}
			(eventHeaps.hostLength())[threadId] = tmpEventHeap.size();
		}
	
	sites.memcpyH2D();
	centers.memcpyH2D();
	headArcsPos.memcpyH2D();
	bounds.memcpyH2D();
	edges.memcpyH2D();
	vertices.memcpyH2D();
	
	
	iterateTimes.memcpyH2D();
	
	eventHeaps.memcpyH2D();
	arcHeaps.memcpyH2D();
	
	
	
}


void copyArc(){
	for(int threadNo = 0; threadNo < blockNum* threadNum; threadNo++ ){
		for(int i = 0; i < (arcs.hostLength())[threadNo]; ++i){
			EllipseArc* tmpArc = &(arcs.atHost(threadNo,i));
			EllipseArcNoPt* tmpArcNoPt = &(arcsNoPt.atHost(threadNo,i));
			if( tmpArcNoPt->leftActiveEdgePos == -1 ){
				tmpArc->leftActiveEdge = NULL;
			}else{
				tmpArc->leftActiveEdge = &(activeEdges.atHost(threadNo,tmpArcNoPt->leftActiveEdgePos));
			}
			if( tmpArcNoPt->rightActiveEdgePos == -1){
				tmpArc->rightActiveEdge = NULL;
			}else{
				tmpArc->rightActiveEdge = &(activeEdges.atHost(threadNo,tmpArcNoPt->rightActiveEdgePos));
			}
		}
	}
}

void copyActiveEdge()
{
	for(int threadNo = 0; threadNo < blockNum * threadNum; threadNo++){
		for(int i = 0; i < (activeEdgesNoPt.hostLength())[threadNo];++i ){
			ActiveEdge* tmpActiveEdge = &(activeEdges.atHost(threadNo,i));
			ActiveEdgeNoPt* tmpActiveEdgeNoPt = &(activeEdgesNoPt.atHost(threadNo,i));
			if( tmpActiveEdgeNoPt->leftArcPos == -1){
				tmpActiveEdge->leftArc = NULL;
			}else{
				tmpActiveEdge->leftArc = &(arcs.atHost(threadNo,tmpActiveEdgeNoPt->leftArcPos));
			}
			if( tmpActiveEdgeNoPt->rightArcPos == -1){
				tmpActiveEdge->rightArc = NULL;
			}else{
				tmpActiveEdge->rightArc = &(arcs.atHost(threadNo,tmpActiveEdgeNoPt->rightArcPos));
			}

		}
	}
}
void copyEvent()
{
	for(int threadNo = 0; threadNo < blockNum * threadNum;threadNo++){
		for(int i = 0; i < (eventHeaps.hostLength())[threadNo];++i){
			Event* tmpEvent = &(eventHeaps.atHost(threadNo,i));
			EventNoPt* tmpEventNoPt = &(eventHeapsNoPt.atHost(threadNo,i));
			if( tmpEventNoPt->arcPos == -1){
				tmpEvent->arc = NULL;
			}else{
				tmpEvent->arc = &(arcs.atHost(threadNo,tmpEventNoPt->arcPos));
			}
			if( tmpEventNoPt->leftActiveEdgePos == -1){
				tmpEvent->leftActiveEdge = NULL;
			}else{
				tmpEvent->leftActiveEdge = &(activeEdges.atHost(threadNo,tmpEventNoPt->leftActiveEdgePos));
			}
			if( tmpEventNoPt->rightActiveEdgePos == -1){
				tmpEvent->rightActiveEdge = NULL;
			}else{
				tmpEvent->rightActiveEdge = &(activeEdges.atHost(threadNo,tmpEventNoPt->rightActiveEdgePos));
			}
		}
	}
}



void voronoi_kernel_c(){
	cuPrintInit();

	
	voronoi_kernel<<<blockNum,threadNum>>>(sites.devicePoint(),sites.deviceLength(),eventHeaps.devicePoint(),eventHeaps.deviceLength(),vertices.devicePoint(),vertices.deviceLength(),centers.devicePoint(),arcHeaps.devicePoint(),arcHeaps.deviceLength(),arcs.devicePoint(),arcs.deviceLength(),activeEdges.devicePoint(),activeEdges.deviceLength(),edges.devicePoint(),edges.deviceLength(),iterateTimes.devicePoint(),headArcsPos.devicePoint(),arcsNoPt.devicePoint(),arcsNoPt.deviceLength(),activeEdgesNoPt.devicePoint(),activeEdgesNoPt.deviceLength(),eventHeapsNoPt.devicePoint(),eventHeapsNoPt.deviceLength() );


	cutilCheckMsg("kernel launch failure");
#ifdef _DEBUG
//	cutilSafeCall( cudaDeviceSynchronize() );
#endif
	
	sites.memcpyD2H();
	centers.memcpyD2H();
	headArcsPos.memcpyD2H();
	bounds.memcpyD2H();
	
	vertices.memcpyD2H();

	arcs.memcpyD2H();
	arcsNoPt.memcpyD2H();
	copyArc();

	edges.memcpyD2H();
	
	activeEdges.memcpyD2H();
	activeEdgesNoPt.memcpyD2H();
	copyActiveEdge();

	eventHeaps.memcpyD2H();
	eventHeapsNoPt.memcpyD2H();
	copyEvent();

	arcHeaps.memcpyD2H();
	
	cudaPrintfDisplay(stdout, true);
	//printf("after %f\n" ,sites.h_point[0].x );
	//printf("length %d\n", sites.h_length[0]);
	
	cudaPrintfEnd();  
	
}

__global__ void voronoi_kernel(float2* d_sitesOrigin ,int* d_sitesOriginLength,Event* d_eventHeap,int* d_eventHeapSize,float2* d_vertexPool,int* d_vertexPoolLength,float2* d_center,ArcHeapEle* d_arcHeap,int* d_arcHeapSize,EllipseArc* d_arcsPool , int* d_arcsPoolLength,ActiveEdge* d_activeEdgePool,int* d_activeEdgePoolLength,int2* d_edgePool,int* d_edgePoolLength,int*  d_iterateTimes,int* d_headArcPos,EllipseArcNoPt* d_arcsNoPtPool,int* d_arcsNoPtPoolLength,ActiveEdgeNoPt* d_activeEdgeNoPtPool,int* d_activeEdgeNoPtPoolLength,EventNoPt* d_eventHeapNoPt,int* d_eventHeapNoPtSize ){
		
#ifdef _KERNEL_INFO_
	cuPrintf("kernel's block :%d\n", blockDim.x);  
	cuPrintf("kernel's thread :%d\n", blockIdx.x);  
	cuPrintf("kernel's thread :%d\n", blockDim.x*blockIdx.x+threadIdx.x);  
#endif

	int threadNo = blockDim.x * blockIdx.x + threadIdx.x;
#ifndef _DELETE_IT_
	int threadNo = 0;
#endif

		Heap <Event>eventHeap;
		Heap <ArcHeapEle> arcHeap;
		eventHeap.initial(&d_eventHeap[threadNo*EVENT_HEAP_SIZE],d_eventHeapSize[threadNo],EVENT_HEAP_SIZE);
		arcHeap.initial(&d_arcHeap[threadNo*ARC_HEAP_SIZE], d_arcHeapSize[threadNo],ARC_HEAP_SIZE);

		Pool <EventNoPt>eventHeapNoPt;
		eventHeapNoPt.initial(&d_eventHeapNoPt[threadNo*EVENT_HEAP_SIZE],d_eventHeapNoPtSize[threadNo],EVENT_HEAP_SIZE);

		Pool<float2> vertexPool;
		vertexPool.initial(&d_vertexPool[threadNo*VERTEX_POOL_SIZE],d_vertexPoolLength[threadNo],VERTEX_POOL_SIZE);
		Pool<EllipseArc>arcsPool;
		arcsPool.initial(&d_arcsPool[threadNo*ARC_POOL_SIZE],d_arcsPoolLength[threadNo],ARC_POOL_SIZE);
		Pool<EllipseArcNoPt>arcsNoPtPool;
		arcsNoPtPool.initial(&d_arcsNoPtPool[threadNo*ARC_POOL_SIZE],d_arcsNoPtPoolLength[threadNo],ARC_POOL_SIZE);

		Pool<ActiveEdge>activeEdgePool;
		activeEdgePool.initial(&d_activeEdgePool[threadNo*ACTIVE_EDGE_POOL_SIZE],d_activeEdgePoolLength[threadNo],ACTIVE_EDGE_POOL_SIZE);
		Pool<ActiveEdgeNoPt>activeEdgeNoPtPool;
		activeEdgeNoPtPool.initial(&d_activeEdgeNoPtPool[threadNo*ACTIVE_EDGE_POOL_SIZE],d_activeEdgeNoPtPoolLength[threadNo],ACTIVE_EDGE_POOL_SIZE);
		Pool <int2> edgePool;
		edgePool.initial(&d_edgePool[threadNo*EDGE_POOL_SIZE],d_edgePoolLength[threadNo],EDGE_POOL_SIZE);

		//pEllipseArc headArc = d_headArc[threadNo];
		int headArcPos = d_headArcPos[threadNo];
		EllipseArc* headArc;
		if( headArcPos == -1 )
			headArc == NULL;
		else
			headArc = &arcsPool[headArcPos];
		

		float2 center = d_center[threadNo];
		float2* d_sites = &d_sitesOrigin[threadNo*SITES_SIZE];
		float sitesLength = d_sitesOriginLength[threadNo];

#ifdef PRINT_LOG
		cuPrintf("eventtop %f\n" , eventHeap.top().distance);
#endif
		//for(int i = 0; i < d_iterateTimes[threadNo];++i){
#ifdef PRINT_LOG
			cuPrintf("\n*********** total i %d\n" , i );
			cuPrintf("\n****************************\ntotal i %d\n" , i);
#endif
			removeOutofdateEvents(eventHeap);
			if(eventHeap.empty())  return;


			//for(int j = 1; j <= d_eventHeapSize;++j){
			//cuPrintf("d_eventheap %f\n",d_eventHeap[j].distance);
			//}
#ifdef PRINT_LOG

			cuPrintf("distance %f \n",eventHeap.top().distance);
#endif

			
			if( headArc == NULL){
			
				float dOfFirstLevel = eventHeap.top().distance;

				Event lastEvent;
				lastEvent = eventHeap.top();

				int cnt = 0;
				while( eventHeap.top().type == NewEllipse &&
					eventHeap.top().distance - dOfFirstLevel < WXN_EPSILON ){
						cnt++;
						EllipseArc* arc;
						
						arcsPool.alloc(arc);
						arc->initial(center,eventHeap.top().site,NULL,NULL);

						if( headArc == NULL ){
							headArc = arc;
						}
						ArcHeapEle tmpArcHeapEle(arc->GetRotateAngle(),arc);
						arcHeap.push(tmpArcHeapEle);
						lastEvent = eventHeap.pop();
				}

				

				eventHeap.push(lastEvent);


			}else if(eventHeap.top().type == NewEllipse){
#ifdef PRINT_LOG
				//cuPrintf("y %f , x %f\n",heapTop<Event>(d_eventHeap,d_eventHeapSize)->site.y - center.y,heapTop<Event>(d_eventHeap,d_eventHeapSize)->site.x - center.x);
				cuPrintf("y %f , x %f\n" , eventHeap.top().site.y - center.y
					,eventHeap.top().site.x - center.x);
#endif
				float angleSite = atan2f(eventHeap.top().site.y - center.y
					,eventHeap.top().site.x - center.x);
				EllipseArc* arc = headArc;
				bool find = false;
				bool fEntireCircle = (headArc->leftActiveEdge==NULL);
				bool fCloseToLeft = false;
				bool fCloseToRight = false;

#ifdef PRINT_LOG

				cuPrintf("-------------------\nangleSite:%f\n",angleSite);
				do{
					//float2 angles = arc->GetAngles(heapTop<Event>(d_eventHeap,d_eventHeapSize)->distance);
					float2 angles = arc->GetAngles(eventHeap.top().distance);
					cuPrintf("angles : (%f , %f)\n",angles.x , angles.y);
				}while(arc->leftActiveEdge != NULL && ( arc = arc->leftActiveEdge->leftArc ) != headArc );
				arc = headArc;
				cuPrintf("----------------------\n");
#endif

				
				int tmpcnt = 0;
				do{
					float2 angles = arc->GetAngles(eventHeap.top().distance);

					float tmpAngle = angleSite;
#ifdef PRINT_LOG

					cuPrintf("tmpAngle %f",tmpAngle);
#endif

					if( tmpAngle < angles.x )
						tmpAngle += 2 * M_PI;
					if( tmpAngle >= angles.x && tmpAngle <= angles.y ){
						find = true;
						if( tmpAngle - angles.x < WXN_EPSILON )
							fCloseToRight = true;
						if( tmpAngle - angles.y > -WXN_EPSILON )
							fCloseToLeft = true;
						break;
					}
					if( arc->leftActiveEdge == NULL )
						arc = NULL;
					else
						arc = arc->leftActiveEdge->leftArc;
				}while( arc != NULL && arc != headArc );

				

				if( !find ){
					cuPrintf("arc not find!");
					return;
				}
				EllipseArc* arcNew;
				arcsPool.alloc(arcNew);
				arcNew->initial(center,eventHeap.top().site,NULL,NULL);
				float2 intersection;



				if(arc->GetP(eventHeap.top().distance) < WXN_EPSILON){
					intersection = center;
				}else{
					intersection = arc->GetPointAtAngleParam(eventHeap.top().distance,angleSite);
				}
				vertexPool.insert(intersection);


				

				int nV = vertexPool.size() - 1;
				if(fEntireCircle){
					//cuPrintf("entire circle");
					
					ActiveEdge* activeEdge1;
					activeEdgePool.alloc(activeEdge1);
					if( !activeEdge1->initial(arc,arcNew,nV) ) return;
					ActiveEdge* activeEdge2;
					activeEdgePool.alloc(activeEdge2);
					if( !activeEdge2->initial(arcNew,arc,nV) ) return;

					arc->leftActiveEdge = activeEdge2;
					arc->rightActiveEdge = activeEdge1;
					arcNew->leftActiveEdge = activeEdge1;
					arcNew->rightActiveEdge = activeEdge2;
					
				}else{
					if(fCloseToLeft){
						EllipseArc* left = arc->leftActiveEdge->leftArc;
						int2 edge = make_int2(left->rightActiveEdge->fixedVertex,nV);
						edgePool.insert(edge);
						activeEdgePool.alloc(left->rightActiveEdge);
						if( !left->rightActiveEdge->initial(left,arcNew,nV) ) return;
						ActiveEdge* activeEdge;
						activeEdgePool.alloc(activeEdge);
						if( !activeEdge->initial(arcNew,arc,nV) ) return;
						arcNew->leftActiveEdge = left->rightActiveEdge;
						arcNew->rightActiveEdge = activeEdge;
						arc->leftActiveEdge = activeEdge;
						Event evt;
						left->GetVanishEvent(evt,vertexPool);
						if( evt.distance < FLT_MAX ){
							eventHeap.push(evt);
						}
						arcNew->GetVanishEvent(evt,vertexPool);

						if( evt.distance < FLT_MAX ){
							eventHeap.push(evt);
						}
						arc->GetVanishEvent(evt,vertexPool);
						if( evt.distance < FLT_MAX ){
							eventHeap.push(evt);
						}
					}else if( fCloseToRight ){

						
						EllipseArc* right = arc->rightActiveEdge->rightArc;
						int2 edge = make_int2(right->leftActiveEdge->fixedVertex,nV);
						edgePool.insert(edge);
						activeEdgePool.alloc(right->leftActiveEdge);
						if( !right->leftActiveEdge->initial(arcNew,right,nV) ) return;
						ActiveEdge* activeEdge;


						activeEdgePool.alloc(activeEdge);


						if( !activeEdge->initial(arc,arcNew,nV) ) return;

						arcNew->rightActiveEdge = right->leftActiveEdge;
						arcNew->leftActiveEdge = activeEdge;
						arc->rightActiveEdge = activeEdge;
						Event evt;
						right->GetVanishEvent(evt,vertexPool);
						if( evt.distance < FLT_MAX ){
							eventHeap.push(evt);
						}
						arcNew->GetVanishEvent(evt,vertexPool);
						if( evt.distance < FLT_MAX ){
							eventHeap.push(evt);
						}
						arc->GetVanishEvent(evt,vertexPool);
						if( evt.distance < FLT_MAX ){
							eventHeap.push(evt);
						}
					}else{
						
						EllipseArc* arc1 = arc;
						ActiveEdge* rightActiveEdge = arc->rightActiveEdge;
						EllipseArc* arc2;
						//if( !allocArc( arc2,d_arcsPool,d_arcsPoolLength) ) return;
						arcsPool.alloc(arc2);
						*arc2 = *arc;
						ActiveEdge* activeEdge1;
						//if( !allocActiveEdge( activeEdge1,d_activeEdgePool,d_activeEdgePoolLength) ) return;
						activeEdgePool.alloc(activeEdge1);
						if( !activeEdge1->initial(arc1,arcNew,nV) ) return ;
						ActiveEdge* activeEdge2;
						//if( !allocActiveEdge( activeEdge2,d_activeEdgePool,d_activeEdgePoolLength) ) return;
						activeEdgePool.alloc(activeEdge2);
						if( !activeEdge2->initial(arcNew,arc2,nV) ) return;
						arc2->leftActiveEdge = activeEdge2;
						arc1->rightActiveEdge = activeEdge1;
						arcNew->leftActiveEdge = activeEdge1;
						arcNew->rightActiveEdge = activeEdge2;
						arc2->rightActiveEdge = rightActiveEdge;
						if(rightActiveEdge != NULL){
							rightActiveEdge->leftArc = arc2;
						}
						Event evt;
						arc1->GetVanishEvent(evt,vertexPool);
						if( evt.distance < FLT_MAX ){
							eventHeap.push(evt);
						}
						arc2->GetVanishEvent(evt,vertexPool);
						if( evt.distance < FLT_MAX ){
							eventHeap.push(evt);
						}
					}
				}

				
			} else if (eventHeap.top().type == EllipseVanish){
				float actualDis = eventHeap.top().distance + TOLERANCE;
				if( eventHeap.top().leftActiveEdge == eventHeap.top().arc->leftActiveEdge &&
					eventHeap.top().rightActiveEdge == eventHeap.top().arc->rightActiveEdge){
						float2 v = eventHeap.top().site;
						vertexPool.insert(v);
						int vIndex = vertexPool.size() - 1;
						headArc = eventHeap.top().leftActiveEdge->leftArc;

						EllipseArc* inactiveArcs;

						EllipseArc* leftArc = eventHeap.top().leftActiveEdge->leftArc;
						EllipseArc* rightArc = eventHeap.top().rightActiveEdge->rightArc;
						edgePool.insert(make_int2(eventHeap.top().arc->leftActiveEdge->fixedVertex,vIndex));
						edgePool.insert(make_int2(eventHeap.top().arc->rightActiveEdge->fixedVertex,vIndex));

						inactiveArcs = eventHeap.top().arc;

						ActiveEdge* newActiveEdge;
						//allocActiveEdge(newActiveEdge,d_activeEdgePool,d_activeEdgePoolLength);
						activeEdgePool.alloc(newActiveEdge);
						if( !newActiveEdge->initial(leftArc,rightArc,vIndex) ) return ;
						leftArc->rightActiveEdge = newActiveEdge;
						rightArc->leftActiveEdge = newActiveEdge;

						inactiveArcs->leftActiveEdge = inactiveArcs->rightActiveEdge = NULL;
						if(inactiveArcs == headArc){
							headArc = rightArc;
						}


						Event evt;
						leftArc->GetVanishEvent(evt,vertexPool);
						if(evt.distance < FLT_MAX){
							eventHeap.push(evt);
						}
						rightArc->GetVanishEvent(evt,vertexPool);
						if(evt.distance <FLT_MAX){
							eventHeap.push(evt);
						}
				}
			}
			
			eventHeap.pop();

		//}


		nxt:
#ifdef _KERNEL_INFO_
		cuPrintf("d_eventHeapsize %d\n" , eventHeap.maxEverSize());
		cuPrintf("arcHeapSize %d\n" , arcHeap.maxEverSize());
		cuPrintf("d_vertexPoolLength %d\n" , d_vertexPoolLength[threadNo]);
		cuPrintf("d_arcsPoolLength %d\n" , d_arcsPoolLength[threadNo]);
		cuPrintf("d_activeEdgePoolLength %d\n" ,d_activeEdgePoolLength[threadNo]);
		cuPrintf("d_edgePoolLength %d\n" , d_edgePoolLength[threadNo]);
#endif
		d_eventHeapSize[threadNo] = eventHeap.size();
		d_arcHeapSize[threadNo] = arcHeap.size();
		d_vertexPoolLength[threadNo] = vertexPool.size();
		d_arcsPoolLength[threadNo] = arcsPool.size();
		d_activeEdgePoolLength[threadNo] = activeEdgePool.size();
		d_edgePoolLength[threadNo] = edgePool.size();
		d_center[threadNo] = center;

		
		d_eventHeapNoPtSize[threadNo] = d_eventHeapSize[threadNo];
		for(int i = 0; i < eventHeap.size();++i){
		 	eventHeapNoPt[i].initial(&eventHeap[i],&arcsPool[0],&activeEdgePool[0]);
		}

		d_arcsNoPtPoolLength[threadNo] = d_arcsPoolLength[threadNo];
		for(int i = 0; i < arcsPool.size();++i){
			(arcsNoPtPool[i]).initial( &arcsPool[i] , &(activeEdgePool[0]) );
		}
		d_activeEdgeNoPtPoolLength[threadNo] = d_activeEdgePoolLength[threadNo];
		for(int i = 0; i < activeEdgePool.size();++i){
			(activeEdgeNoPtPool[i]).initial(&activeEdgePool[i],&(arcsPool[0]));
		}

		if( headArc == NULL ){
			d_headArcPos[threadNo] = -1;
		}else{

			int i;
			for(i = 0; i < arcsPool.size();++i){
				if( headArc == &arcsPool[i] ){
					//cuPrintf("head %d\n",i);
					break;
				}
			}

			d_headArcPos[threadNo] = i;
		}

}

__device__ __host__ void removeOutofdateEvents(Heap<Event>& eventHeap)
{
	while( !eventHeap.empty() ){
		if( eventHeap.top().distance > FLT_MAX){
			eventHeap.pop();
		}else if( eventHeap.top().type == NewEllipse ){
			break;
		}
		if( eventHeap.top().leftActiveEdge == eventHeap.top().arc->leftActiveEdge &&
			eventHeap.top().rightActiveEdge == eventHeap.top().arc->rightActiveEdge ){
				break;
		}else{
			eventHeap.pop();
		}
	}
}

__device__ __host__ Event GetNextEvent(Heap<Event>& eventHeap)
{
	Event evt;
	evt.distance = FLT_MAX;
	evt.arc = NULL;
	evt.leftActiveEdge = evt.rightActiveEdge = NULL;
	evt.type = -1;
	removeOutofdateEvents(eventHeap);
	if(eventHeap.empty()) return evt;
	evt = eventHeap.top();
	return evt;
}

#endif